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the  theoretical  formulations  and  numerical  simulations  of  electromagnetic  wave  interaction 
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with  adhesive  scatterers,  and  their  applications  in  atmospheric  aerosols  and  random  media 
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Development  of  a  fractal-based  aggregated  model  and  a  rigorous  T-matrix  approach  for 
studying  the  light  scattering  and  absorption  by  agglomerated  particles.  We  have  simulated 
the  structure  of  random  clusters,  determine  the  fractal  dimension,  and  investigated  their 
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optical  properties,  (ii)  Development  of  a  sticky  particle  model  for  simulating  a  random 
medium  consisting  of  aggregated  particles.  We  have  studied  the  propagation  and  scattering 
of  waves  in  such  random  media  and  performed  computer  simulations  to  determine  the 
cluster  size  distribution  and  the  fractal  dimension  for  agglomerates.  In  the  following,  we 
describe  these  major  accomplishments. 
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1.  STATEMENT  OF  RESEARCH  PROBLEM 


* 


The  aerosols  in  atmosphere  are  often  random  cluster  of  small  primary  particles  with 
ramified  appearance.  They  are  considered  to  be  fractal  aggregates  with  a  noninteger  di¬ 
mensions  and  they  have  the  important  property  of  invariance  under  scale  transformation. 
This  essential  fractal  morphology  has  already  been  demonstrated  for  soot  and  smoke.  In 
this  research  we  investigate  the  scattering  and  absorption  properties  of  aggregated  aerosol 
particles  based  on  a  fractal  cluster  model.  A  rigorous  treatment  of  the  scattering  of  elec¬ 
tromagnetic  radiations  is  pursued  to  take  into  account  the  features  of  inter-connected 
primary  particles.  Specific  research  efforts  include  using  the  T-matrix  method  together 
with  the  translation  addition  theorem  for  vector  spherical  waves  to  solve  the  scattering 
and  absorption  by  connected  spherules  with  clustering  and  branched  structures;  computer 
simulating  the  configurations  for  random  aggregates,  solving  the  scattering  and  absorp¬ 
tion  cross  sections  for  random  agglomerates;  calculating  the  average  optical  properties  of 
random  clusters  over  configurations;  and  validating  the  theoretical  model  by  comparison 
with  experimental  measurements. 


2.  SUMMARY  OF  RESULTS 
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The  research  performed  during  the  period  of  July  1991  -  September  1993  was  focused  on 
the  theoretical  formulations  and  numerical  simulations  of  electromagnetic  wave  interaction 
with  random  fractal  clusters  and  electromagnetic  wave  propagation  and  scattering  in  media 
with  adhesive  scatterers,  and  their  applications  in  atmospheric  aerosols  and  random  media 
consisting  of  aggregated  particles.  The  major  accomplishments  in  this  period  include:  (i) 
Development  of  a  fractal-based  aggregated  model  and  a  rigorous  T-matrix  approach  for 
studying  the  light  scattering  and  absorption  by  agglomerated  particles.  We  have  simulated 
the  structure  of  random  clusters,  determine  the  fractal  dimension,  and  investigated  their 
optical  properties,  (ii)  Development  of  a  sticky  particle  model  for  simulating  a  random 
medium  consisting  of  aggregated  particles.  We  have  studied  the  propagation  and  scattering 
of  waves  in  such  random  media  and  performed  computer  simulations  to  determine  the 
cluster  size  distribution  and  the  fractal  dimension  for  agglomerates.  In  the  following,  we 
describe  these  major  accomplishments. 

2.1  Light  Scattering  and  Absorption  by  Agglomerated  Particles 

The  aerosols  in  atmosphere  are  often  random  clusters  of  small  primary  particles  with 
ramified  and  connected  appearance.  They  are  considered  to  be  fractal  aggregates  and 
the  essential  fractal  morphology  has  already  been  demonstrated  for  soot  and  smoke  [1-4]. 
The  aggregation  process  of  small  primary  particles  has  been  simulated  by  applying  the 
cluster-cluster  aggregation  (CCA)  algorithm  [5-10].  In  this  model,  the  primary  particles 
undergo  random  walk  simultaneously  to  form  small  clusters  initially,  and  subsequent  col¬ 
lisions  will  lead  to  their  aggregation.  The  fractal  dimensions  of  resultant  agglomerates  are 
determined  by  finding  out  the  power  form  relationship  between  the  number  of  primary 
particles  and  the  cluster  gyration  radius  [9-10].  To  study  the  scattering  and  absorption 
of  electromagnetic  waves  by  such  fractal  aggregates,  a  T-matrix  approach  together  with 
the  translation  addition  theorem  for  vector  spherical  waves  is  employed  to  solve  the  ag¬ 
gregate  scattering  problem  [11-12].  Numerical  results  for  the  phase  function,  scattering 
and  absorption  cross  sections  are  illustrated  as  a  function  of  observation  angles,  cluster 
structures,  particle  sizes  and  number  of  primary  particles.  The  model  predicted  that  the 
effects  of  multiple  scattering  may  cause  scattering  and  absorption  cross  sections  of  an 
aggregate  that  will  initially  increase  with  its  size  (number  of  primary  particles  or  size  of 
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monomer),  then  saturate  and  decrease  when  its  size  further  increase.  We  also  made  the 
comparison  with  the  measurements  of  extinction  cross  sections  on  soot  particles  over  the 
wavelength  region  from  millimeter  wave  to  infrared  [13].  It  is  shown  that  the  extinction 
decreases  linearly  with  increasing  wavelength,  however,  the  slope  is  not  in  agreement  with 
the  measurements. 

2.1.1  Multiple  Scattering  Formulation 

Consider  a  plane  wave  incident  on  an  aggregate  in  the  direction  of  fc,-,  =  sin0j  cos  4>ii+ 

sin  sin  <f>iy  4-  cos  6{Z.  The  aggregate  is  modeled  as  a  cluster  of  N  connected  identical 
spheres  centered  at  ¥if  r2,  . . .,  r/y,  each  of  which  has  radius  a  and  refractive  index  /*,. 
The  background  medium  has  the  refractive  index  p.  In  a  cluster  of  N  particles,  the 
total  amount  of  wave  excitation  exerted  on  the  particle  t  is  the  sum  of  incident  field  and 
scattered  fields  from  the  other  particles, 


where  and  o*"*  are  column  vectors  that  represent  the  exciting  field  of  the  scatterer  l 
and  the  incident  field  respectively,  T  is  the  T  matrix  that  describes  scattering  from  the 
scatterer  j,  W(rtrj)  is  the  transformation  matrix  used  to  convert  spherical  waves  centered 
at  Tj  to  spherical  waves  centered  at  r>  (see  Appendix  A),  and  k  is  the  wave  number  of  the 
background  medium  [12]. 

The  multiply  scattered  field  5*^  of  particle  t  is  given  as 

3 (2) 

where  o is  the  column  vector  that  represent  the  scattered  field  from  the  particle  i.  The 
total  scattered  field  at  a  field  point  ¥  is  given  by  the  summation  of  scattered  waves  from 
all  particles 

=  E£  + a%Mw^(kmj}  (3) 

4=1  "».n 

where  Mmn  and  7Tmn  are  vector  spherical  waves  (see  Appendix  B)  [12],  and  and 

amn  ^  are  expansion  coefficients  for  the  scattered  field  form  particle  l.  The  far  scattered 
field  in  the  direction  &#,  fc,  =  sin0#  cos  <£,x+sin0,  sin^,y+cos  0tz,  can  be  expressed  as 

F  =  j^(v.Ev.  +  h.Eh.)  (4) 
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t 


with  v,  =  cos9,  cos<f>,x+cos8,  sin  ^,y— sin  6,z  and  ht  =  —  s\n<t>,x+cos<f>,y  being  the  vertical 
and  horizontal  polarization  vectors  respectively  for  the  scattered  wave.  From  (3)  and  (4), 
Ev,  and  Eh,  can  be  expressed  as 


E„  =  £  4.)  +  «•«£?*  •  (5) 

m,n 

Eh.  =  £w‘-'  +  (6) 

m,r» 

where  J5mn  and  Cmn  are  vector  spherical  harmonics  [12],  7mn  are  coefficients  as  defined  in 
Appendix  B,  and 


N 


.(M) 


=  £ 


(7) 


l=\ 

N 


a'(N)  _  y'  *N)(t) 

ttmn  /  ^  “  amn 

*=1 


(8) 


Instead  of  solving  the  exciting  field  coefficients,  we  can  combine  equations  (1)  and  (2) 
to  form  equations  for  the  scattered  field  coefficients 

<?(*)  =  £  •  W(kJ7f])  •  o*<»  +  ¥t)-aine  (9) 

Thus  the  multiple  scattering  solution  consists  of  solving  (9)  numerically  for  a’W  and  then 
using  (3)  to  calculate  the  scattered  field  2?*(F). 

Equation  (9)  is  solved  by  using  an  iterative  scheme  [14],  The  result  for  the  (^  +  1) 
iteration  is 

o*^1)  =  eiS<  ?*  ^  ■  jr*  +  £  •  w(k¥FJ)  •  a*0*00  (io) 

ii* 

where  v  denotes  the  i/th  iterated  solution.  The  initial  solution  of  a*^1*  is  just  the  first  term 
on  the  right-hand  side  of  (10).  The  physical  correspondence  is  the  first-order  scattering, 
the  second-order  scattering,  and  so  on. 

The  bistatic  scattering  cross  section  is  defined  by 

*«.(*.,  =  (11) 

\Ea  I 

where  a  and  0  represent  the  polarization  states  of  the  incident  and  scattered  wave  respec¬ 
tively,  and  a  and  0  can  be  v  or  h. 
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2.1.2  Generation  of  Random  Aggregates 


Aggregation  is  an  irreversible  growth  process  in  which  initially  dispersed  primary  parti¬ 
cles  collide  and  stick  together  to  build  the  clumping  structure.  For  example,  soot  particles 
are  formed  by  random  collision  of  small  carbon  spheres  in  the  turbulent  region  of  a  flame. 
Growth  models  such  as  diffusion-limited  aggregation  (DLA)  [5]  and  cluster-cluster  aggre¬ 
gation  (CCA)  (6]  have  been  applied  to  obtain  information  about  the  growth  processes  and 
predict  the  fractal  structure  of  aggregates.  In  the  DLA  model  primary  particles  are  added 
to  a  growing  cluster  successively.  Each  new  particle  moves  in  a  random  manner,  its  path 
being  chosen  by  the  Monte  Carlo  method.  The  CCA  model  is  a  variation  of  the  DLA 
model,  in  which  many  particles  undergoing  random  walk  simultaneously  to  form  small 
clusters  initially,  and  subsequent  collisions  will  lead  to  their  aggregation.  However,  many 
studies  of  soot  particles  have  shown  that  they  are  formed  by  cluster-cluster  aggregation 
[5-7]  with  an  open  structure  and  the  fractal  dimensions  are  in  the  range  of  1.7  —  1.9. 

The  mechanism  used  here  for  generating  fractal  clusters  is  the  hierarchical  cluster- 
cluster  aggregation  algorithm  [8].  In  this  model  collisions  are  restricted  to  clusters  of 
the  same  size  during  the  simulation,  where  dimers  collide  and  stick  to  form  tetramers, 
tetramers  collide  and  stick  to  make  octamers,  etc.  This  simulation  process  is  more  speedy 
than  the  CCA  model  and  produces  clusters  of  very  similar  fractal  dimensions  [8].  Examples 
of  3-dimensional  aggregates  obtained  by  this  process  are  illustrated  in  Figure  1. 

For  fractal  clusters,  despite  their  random  appearance,  there  exists  a  power  form  rela¬ 
tionship  between  the  number  of  primary  spheres  N  and  the  radius  of  gyration  Rg 

N  =  cR°'  (12) 


where  Df  is  the  fractal  dimension  and  c  is  a  proportionality  constant.  Assuming  all 
primary  spheres  to  have  equal  mass,  the  radius  of  gyration  of  an  aggregate  is  given  by 


R\ 


(13) 


where  is  the  distance  from  the  monomer  *  to  the  center  of  aggregate.  The  fractal 
dimension  D/w  1.74  is  deduced  from  the  ln(iV)-ln(iZ9)  data  curve  fitting  in  Figure  2. 


2.1.3  Simulation  Results 


In  this  section  we  simulate  the  scattering  of  optical  waves  by  soot  aggregates  that  are 
generated  from  the  hierarchical  CCA  model  with  ka  =  0.05  and  ka  =  0.25.  A  typical 
refractive  index  for  soot  particle  /z#  =  1.75+ *0.5  is  considered  in  the  simulation. 


Assuming  k <  =  £,  we  calculate  the  bistatic  scattering  coefficients  as  a  function  of  scat¬ 
tering  angles.  For  Figures  3(a)  and  3(b),  we  consider  the  scattering  by  a  cluster  of  N  =  32 
primary  spheres  with  ka  =  0.05  and  ka  =  0.25,  respectively.  We  plot  bistatic  scattering 
coefficients  as  a  function  of  <j>t  with  9,  =  90°.  The  coefficients  show  that  a  maximum  scat¬ 
tering  is  in  the  forward  direction  and  the  result  for  vv  is  larger  than  that  for  hh.  There  is 
a  minimum  for  hh  around  <f>,  =  90°  because  the  induced  dipoles  in  small  primary  spheres 
are  pointing  in  the  same  direction  as  kt.  In  Figures  4(a)  and  4(b),  we  illustrate  the  case 
for  a  larger  cluster  with  N  =  128  monomers.  Figure  4(a)  shows  that  for  a  larger  cluster  we 
have  a  much  stronger  forward  scattering  than  that  in  Figure  3.  The  minima  around  90° 
for  the  hh  are  also  due  to  the  induced  dipoles’  orientations. 


In  Figure  5,  we  plot  the  ratio  of 


NaP 


as  a  function  of  the  number  of  primary  particles 


in  a  cluster,  a W  is  the  scattering  cross  section  for  an  agglomerate  with  N  number  of 
monomers,  a J1)  is  the  scattering  cross  section  for  a  primary  particle.  The  results  have  been 
averaged  over  30  realizations.  It  shows  that  the  scattering  increases  with  N  when  clusters 
are  small  (small  N  or  small  ka).  However,  the  multiple  scattering  becomes  important  as 
clusters  grow  larger  which  is  indicated  by  the  saturation  of  scattering  cross  section.  The 


ratio  of 


N<ril) 


is  shown  in  Figure  6,  where  the  subscript  a  is  used  to  denote  the  absorption 


cross  section.  As  expected,  the  absorption  increases  with  N  when  multiple  scattering  is 
a  minor  effect.  However,  when  the  agglomerates  become  larger,  the  multiple  scattering 
effects  will  block  the  absorption  for  the  inner  particles  within  the  aggregates.  In  Figure  7, 
we  illustrate  the  results  for  extinction  cross  section  which  is  equal  to  the  sum  of  scattering 
and  absorption  cross  sections  <re  =  cr,+<Ta. 


Figure  8  shows  the  wavelength  dependence  of  extinction  cross  section  per  constituent 
particle.  The  index  of  refraction  for  carbon  as  a  function  of  frequency  is  obtained  from  [15]. 
It  is  shown  that  the  extinction  decreases  linearly  with  increasing  wavelength,  but  the  slope 
is  not  in  agreement  with  the  experimental  measurements  [13].  However,  in  the  above  work 
we  have  assumed  only  dipole  interaction  between  primary  particles  and  they  touching  only 
at  one  point  on  the  surface.  This  assumption  has  been  justified  for  dielectric  particles, 
because  in  this  cue  the  displacement  current  dominates  over  the  conduction  current.  In 
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soot  particles,  the  carbon  monomers  may  overlap  by  nucleation  and  may  have  very  high 
conductivity.  For  these  cases,  we  may  need  to  include  more  multipole  moments  in  the 
T-matrix  calculation  to  take  into  account  the  higher  order  interaction  between  primary 
particles.  Also,  the  sphere  model  may  not  be  an  accurate  one  for  highly  absorptive  cases. 
This  is  due  to  the  difficulty  of  field  penetration  into  the  sphere. 

2.2  Scattering  of  Waves  in  Media  with  Adhesive  Particles 

Some  terrain  media  (e.g.  snow,  ice)  may  consists  of  aggregated  particles,  also  they 
can  be  characterized  as  dense  media.  In  a  dense  medium,  the  inhomogeneities  occupy 
an  appreciable  volume  fraction,  the  classical  assumption  of  independent  scattering  is  not 
valid,  we  have  to  consider  the  effects  of  correlated  scattering  and  mutual  coherent  wave 
interactions  among  scatterers.  However,  the  aerosol  can  be  characterized  as  locally  dense 
medium  because  of  the  primary  particles  being  aggregated  together  within  a  cluster.  In  or¬ 
der  to  account  for  the  clustering  of  particles,  we  have  applied  the  sticky  hard  sphere  (SHS) 
model  [16]  by  introducing  a  surface  adhesive  force  between  two  particles.  The  analytic 
wave  approaches,  quasicrystalline  approximation  (QCA)  and  quasicrystalline  approxima¬ 
tion  with  coherent  potential  (QCA-CP)  [12],  which  take  into  account  the  mutual  coherent 
wave  interactions,  were  applied  to  derive  closed-form  equations  for  effective  propagation 
constants  in  random  media  with  sticky  particles.  The  model  produced  results  that  explain 
the  observed  features  of  the  laboratory  experiment  [17].  The  experiment  involved  light 
scattering  from  colloidal  silica  particles  suspended  in  a  solvent  background.  The  compari¬ 
son  with  theory  shows  that  the  experimental  data  can  be  explained  by  assuming  the  silica 
particles  have  some  adhesion,  but  cannot  explained  without  this  assumption.  The  sticky 
particle  pair  distribution  function  has  been  used  by  other  researchers  to  explain  x-ray 
scattering  from  aggregated  gold  colloids  [18].  We  also  performed  Monte  Carlo  simulations 
to  generate  random  media  with  aggregated  particles.  The  computer  simulated  random 
aggregated  media  has  been  shown  to  have  fractal  structure,  and  the  fractal  dimension 
depends  on  the  concentration  and  stickiness  of  particles. 

2.2.1  Pair  Distribution  Function  for  Media  with  Sticky  Particles 

We  shall  consider  a  system  consisting  of  non-interpenetrable,  spherical  particles  of 
diameter  d  with  a  non-zero  surface  adhesive  force,  i.e.  sticky  hard  spheres  (SHS).  In 
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this  model,  the  interaction  between  two  particles  is  of  very  short  range,  the  interparticle 
potential  u(r)  is  described  by  [16]: 

{oo  for  0  <  r  <  s 

In  12t^~*)  for  3  <  r  <  d  (14) 

0  for  r  >  d 

where  the  limit  is  taken  in  such  a  way  that  the  relation 

lim(d  -  a)e~u(r)  =  ^  <  °°  (15) 

holds  with  d  and  r  being  held  fixed.  The  parameter  r  in  (14)  is  dimensionless,  and  its 
inverse  is  a  measure  of  the  attraction  or  stickiness  between  particles.  The  case  of  r-1  =  oo 
corresponds  to  infinite  stickiness  and  r-1  =0  corresponds  to  non-sticky  particles. 

The  pair  distribution  function  g(r)  measures  the  probability  of  finding  a  particle  at  a 
point  r  given  a  particle  at  r  =  0.  The  total  correlation  function  h(r)  between  a  pair  of 
particles  is  defined  as 

Mr)  =  s(r)  - 1  (16) 

The  direct  correlation  function  c(r)  between  a  pair  of  particles,  which  is  short-ranged,  is 
related  to  h(r)  by  means  of  the  Ornstein-Zernike  relation: 

h(r)  =  c(r)  +  n  J  dr*  c(r')  h(\¥  -  r'l)  (17) 

where  n  is  the  number  of  particles  per  unit  volume.  Equation  (17)  indicates  that  the  total 
influence  of  a  particle  on  another  particle  in  the  presence  of  remaining  particles  can  be 
decomposed  into  a  sum  of  the  direct  effect  and  the  indirect  effect  through  other  particles. 

Under  the  Percus-Yevick  (PY)  approximation  [19],  correlation  functions  c(r)  and  h(r) 
are  approximated  by 

c(r)  =  0  for  r  >  d  (18) 

h{r)  =  —  1  +  6(r  —  d)  for  0  <  r  <  d  (19) 

respectively,  with  t  a  dimensionless  parameter  to  be  determined  later.  The  parameter  t 
tends  to  zero  in  the  limit  r_1  =  0.  The  PY  approximation  of  the  pair  distribution  function 
^(r)  for  the  sticky  spherical  particles  can  be  solved  analytically  using  the  factorization 
method  of  Baxter  [20]. 
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The  Ornstein-Zernike  relationship  (17)  can  be  Fourier  transformed  to  obtain  a  conve¬ 
nient  algebraic  equation: 

1  -nC(p)  =  {l  +nJ7(p)}“1  (20) 

where  H(p)  and  C(p)  are  the  three  dimensional  Fourier  transform  of  h(f)  and  c(r).  Ac¬ 
cording  to  the  Wiener-Hopf  technique  due  to  Baxter  [20],  the  left-hand  side  of  (20)  can 
be  factorized  into  the  form 

l-nC(p)  =  Q(p)Q(-p)  (21) 

where  Q(p)  is  defined  by 

Q(p)  =  l-2nn  [  dr  eiprQ(r)  (22) 

Jo 

and  Q(r)  is  a  real  function,  Q(r)  =  0  for  r  >  d.  For  sticky  particles,  in  the  range  0  <  r  <  d, 
Q(r )  has  a  closed-form  expression: 


where 


Q(r)  =  A?—  +  Br  +  D 
£ 

(23) 

a  1  +  2/  —  /i 

A=  a-/)* 

(24) 

(-3/  +  ^<i 

2(1  -  /)’ 

(25) 

o=-4-b“+ % 

(26) 

/  =  -nd3 

(27) 

n  =  t/(l  -  /) 

(28) 

For  a  given  volume  fraction  /  and  stickiness  parameter  r,  the  parameter  t  is  determined 
by  the  quadratic  equation  [16,21) 


f  1 2  ,T  ,  /ln  1  +  //2  0 

12 1  "( 


(29) 


Moreover,  a  further  condition  to  determine  the  solution  of  t  is  that  Q(0)  must  be  positive, 
or  [16,21] 

fi  <  1  +  2/  (30) 


The  procedure  for  calculating  the  pair  function  is  as  follows.  Given  the  particle  diame¬ 
ter  d,  particle  concentration  /,  and  particle  stickiness  r,  the  parameter  t  is  first  determined 
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from  (29)-(30).  By  using  the  parameter  t  and  equations  (23)-(28),  Q(r)  for  0<r<d  can 
be  computed.  Then  the  Fourier  transform  Q(p)  is  calculated  by  using  (22),  and  C(p)  is 
evaluated  by  using  (21).  Next,  H(p )  is  computed  by  solving  equation  (20).  Then  k(r)  is 
obtained  by  taking  an  inverse  Fourier  transform  from  H(p).  Finally,  the  pair  distribution 
function  g(r)  can  be  calculated  from  (16). 

In  Figure  9,  the  pair  distribution  functions  are  shown  for  systems  with  sticky  particles 
of  identical  size,  r  =  0.2  and  r  =  0.5,  and  compared  with  the  non-sticky  case  for  /  =  0.3.  The 
major  features  for  the  pair  distribution  functions  of  sticky  particles  are  the  occurrence  of 
discontinuities  at  the  particle  separation  equals  one  and  two  diameters.  The  discontinuity 
at  r  —  d  is  because  the  particles  can  not  interpenetrate  each  other.  The  height  of  the  peak 
at  r  =  d  grows  rapidly  with  the  increasing  of  particle  stickiness  while  the  width  of  this 
peak  reduces,  which  shows  a  stronger  connectedness  between  particles  as  particles  getting 
more  sticky,  and  other  particles  are  more  likely  to  be  excluded  from  the  region  d<r<  2d. 
Therefore,  in  a  dense  medium  with  sticky  particles,  the  sticky  particles  tend  to  aggregate 
together.  The  discontinuity  at  r  —  2d  arises  from  the  fact  that  for  r  >  2d,  physically,  when 
the  separation  between  two  particles  is  larger  than  twice  the  diameter,  the  probability  of 
these  two  particles  bound  or  connect  to  a  third  particle  drops  to  zero.  In  Figure  10,  the 
pair  distribution  functions  are  plotted  for  /  =  0.2  and  /  =  0.4,  with  stickiness  r  =  0.2. 
For  higher  concentration  of  sticky  particles,  the  pair  distribution  function  displays  more 
fluctuations  just  as  the  case  of  non-sticky  particles. 

2.2.2  Monte  Carlo  Simulation  of  Media  with  Sticky  Particles 

Basically,  the  Monte  Carlo  technique  used  in  simulating  the  clustering  of  particles  is  to 
deposit  particles  one  by  one  into  a  box  which  contains  all  particles  subject  to  the  condition 
that  no  overlap  between  particles  is  adlowed  [22].  However,  the  essence  in  the  simulation 
process  for  the  case  of  sticky  particles  is  that  we  have  to  calculate,  before  depositing  each 
trial  particle,  the  probabilities  that  determine  the  bound  or  unbound  states  for  the  triad 
particle.  These  probabilities  depend  on  the  stickiness  of  particles,  the  size  of  particle,  and 
the  number  of  particles  existing  in  the  box  [23-24].  After  selecting  a  particular  bound  or 
unbound  state  for  the  trial  particle,  we  will  determine  its  possible  position  in  accordance 
with  the  selected  bound  state  according  to  a  uniform  probability.  Then  we  test  for  overlap 
and  reject  the  triad  position  when  a  overlap  is  detected. 
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Volume  Fraction 

0.2 

0.2 

0.2 

0.3 

0.3 

0.3 

Stickiness  Parameter 

0.2 

0.5 

1 

0.2 

0.5 

1 

Fractal  Dimension 

1.72 

1.74 

1.59 

1.90 

1.76 

1.60 

Constant  c 

1.86 

1.64 

1.73 

1.75 

1.65 

1.84 

Table  1:  Numerical  Results  for  Df  and  c 


In  Figure  11,  the  Monte  Carlo  simulation  of  the  pair  distribution  function  is  shown 
for  a  system  with  N  =  512  sticky  particles  for  the  case  of  /  =  0.3,  and  r  =  0.2.  It  shows 
that  the  Monte  Carlo  simulation  also  indicates  the  discontinuity  in  the  pair  function  when 
the  distance  equals  two  diameters.  Generally,  the  result  shows  a  good  agreement  with  the 
Percus-Yevick  approximation.  A  typical  realization  of  the  simulated  medium  with  sticky 
particles  is  shown  in  Figure  12.  Clearly,  we  can  see  that  particles  form  aggregates  and  the 
structure  of  the  cluster  is  quite  open.  The  fractal  dimensions  of  these  agglomerates  can  be 
calculated  from  the  power-law  relationship  between  the  radius  of  gyration  of  a  cluster  and 
the  number  of  particles  in  a  cluster  as  discussed  in  Section  2.1.2.  The  fractal  dimension 
depends  on  the  concentration  of  particles  and  the  stickiness  of  particles.  Based  on  our 
simulations,  for  512  particles  and  30  realizations,  the  computed  fractal  dimension  is  less 
than  2  and  larger  than  1.7  as  illustrated  shown  in  Figure  13  and  Table  1. 

2.2.3  Computation  of  Effective  Propagation  Constant 

When  the  scatterers  in  a  medium  are  close  to  each  other,  it  is  important  to  consider 
the  mutual  coherent  wave  interactions  between  scatterers  and  the  effects  of  correlated 
scattering.  For  many  years  we  have  used  analytic  wave  theory  and  approximations,  such  as 
quasi-crystalline  approximation  (QCA)  and  quasi-crystalline  approximation  with  coherent 
potential  (QCA-CP),  to  study  the  electromagnetic  wave  propagation  and  scattering  in 
uniformly  dense  media  [12,25-27].  In  these  two  approaches,  the  mutual  coherent  wave 
interactions  among  particles  are  weighted  by  the  pair  distribution  functions.  We  have  also 
applied  Monte  Carlo  method  to  directly  solve  the  Foldy-Lax  multiple  scattering  equations 
for  media  with  densely  distributed  particles.  The  results  of  Monte  Carlo  simulations 
agree  very  well  with  those  of  QCA  and  QCA-CP  [14],  which  is  also  confirmed  by  using 
the  recursive  T-matrix  numerical  algorithm  [28].  Also,  controlled  laboratory  experiments 
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have  verified  the  QCA  and  QCA-CP  theories  [29].  The  important  advantage  of  using  these 
two  analytic  approximations  is  that,  since  they  start  with  Maxwell’s  equations,  the  effects 
of  correlated  scattering  from  different  scatterers,  the  mutual  coherent  wave  interactions 
between  scatterers,  and  spatial  correlation  among  particles  are  all  included  in  their  analytic 
formulations.  Moreover,  to  solve  the  scattering  problem  in  a  random  medium  using  these 
two  analytic  approaches,  we  only  need  to  solve  the  integral  equation  of  the  average  field 
once,  which  is  computationally  much  more  efficient  than  the  Monte  Carlo  approach  which 
needs  to  solve  integral  equations  for  many  realizations. 


In  this  section,  the  effective  propagation  constants  for  QCA  and  QCA-CP  are  given, 
and  the  solution  with  the  sticky  particle  pair  function  are  illustrated.  The  dispersion 
relationships  for  QCA  and  QCA-CP  are  given  as  [12,25-27]: 


K2  =  fc2  + 


K2  =  fc2  + 


/(*?  -  *8)  / 
1  +  %£(W)  I 

M  ~ fc2)  f 


1  + 1 


2(k)  -  fc2)fco3(l  -  /)4 


1  +  t 


9  [l  +  -  /)]  [1  +  2/  -  </(l  -  /)]■ 

2(fc,2  -  k2)Ka\l  -  / )4 


) 


(31) 


(32) 


9  [l  +  ^#(1  -  /)]  (1  +  2/  -  if {l  -  f)Y  ) 

respectively,  where  K  is  the  effective  propagation  constant,  /  is  the  fractional  volume  of 
particle,  fc,  is  the  wave  number  inside  the  particle,  fc  is  the  wave  number  of  background 
medium,  a  is  the  radius  of  particle,  and  t  is  the  parameter  in  term  of  stickiness  r  determined 
form  equation  (29). 


The  attenuation  rates  (2Im( if ))  for  various  values  of  stickiness  parameter  r  as  a  func¬ 
tion  of  fractional  volume  are  shown  in  Figure  14.  It  can  be  seen  that  the  sticky  model 
predicts  larger  attenuation  for  all  fractional  volumes  than  the  non-sticky  model.  We  also 
note  that  the  maximum  attenuation  shifts  with  stickiness  which  indicates  the  scattering 
by  effectively  larger  particles  formed  by  the  clustering.  At  a  value  of  r  =  0.2,  the  predicted 
scattering  may  be  greater  than  the  independent  scattering. 


The  concept  of  sticky  particles  has  been  explored  by  Penders  and  Vrij  in  an  inter¬ 
pretation  of  turbidity  (attenuation  of  a  light  beam  by  scattering  when  passing  through 
a  sample)  studies  on  colloidal  silica  particles  [17].  Details  of  the  experimental  procedure 
sue  given  in  Jansen  et.  at.  [30].  In  these  experiments,  results  from  light  scattering  off  of 
high  concentrations  of  silica  particles  is  explained  with  the  sticky  particle  pair  function 
described  in  section  2.2.1. 
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Particle 

SE2 

SJ9 

SJ4 

Solvent 

toluene 

toluene 

benzene 

Refractive  index  of  particle 

1.437 

1.433 

1.425 

Refractive  index  of  solvent 

1.494 

Radius  (nm) 

47 

34.5 

23 

Stickiness  Parameter 

1.0 

0.9 

1.5 

Table  2:  Parameter  Values 


The  three  stearylsilica  particle  dispersions  used  by  Jansen  et.  al.  in  his  turbidity 
study  were  modeled  using  the  multiple  scattering  theory  with  adhesive  hard  spheres.  The 
parameters  are  given  in  Table  2,  where  the  values  for  r  and  a  are  those  for  best  agreement 
with  QCA.  The  He-Ne  laser  had  a  free  space  wavelength  of  632  nm,  and  the  fractional 
volume  of  the  silica  scatterers  was  varied  from  0  —  40%. 

Figure  15  shows  the  turbidity  as  a  function  of  fractional  volume  for  the  three  particle 
species,  with  the  stickiness  given  by  the  values  of  t  in  Table  2.  The  results  of  QCA  agree 
quite  well  with  the  experimental  data.  However,  without  the  assumption  of  adhesion,  the 
turbidity  behavior  cannot  be  explained. 
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Appendix  A 


The  translational  matrix  <r(kf)  for  vector  spherical  waves  in  equation  (9)  is  expressed 
in  the  form  of  [12] 


w.  l(kf)  %kf) 

V  '  B(kf)  A(kr) 


(A.l) 


where  the  matrices  A{kf)  and  B(kr)  are  of  dimension  XmoxX-^mox-  ^ma*  =  Afmox(iVTnai+2) 
and  Nmax  is  the  maximum  order  kept  in  the  summation.  The  expressions  of  their  elements, 
AIU/mn(kr)  and  Bpvrnn{kf ),  are  given  as  follows. 


Anvmn (kr)  =  — (-iyi52a(ro»nl 

7m*  p 

x  hp(kr)Ypm-“(8,  <f>) 

BMim„(kr)  =  — (-l)M+1^a(m,n|-/i,*/|p,p-l)6(n,*/,p) 
7m*  p 

x  hp(kr)Ypm-»(9,<l>) 


where 


a(m,nMp)  =  (-1)^  +  1)  f<"  +  ^  +  ~  ~ 

{n  —  m)\(v  —  ft)l(p  +  m  +  fty. 

x  (  H  U  P  Wn,/P\ 

\  m  ft  -( m  +  n )  )  \  0  0  0  J 


n  v  p 
m  ft  — (m  +  ft) 


a{m,n\fi,v\piq)  =  (-l)m+,1(2p  +  1) 


(n  -f  m)!(i/  +  pt)\(p  ~m-  fi)\  1/2 
(n  -  m)!( v  -  pt)l(p  +  m+  pt)\ 


m  ft  — (m  +  ft) 


)(:::) 


(A.2) 


(A-3) 


(A.4) 


(A.5) 


f-n+p 

a(n,i/,p)=  [  2t'(*/  +  1)(2*^  +  !)  +  (*'  4-  l)(n  +  */ -  p)(n  +  p  -  t/ +  1) 


-  i/(n  +  v  +  p  4-  2)(v  +p  -  n  +  1)]1/J 

6(n,  v,p)  =  ~  2u(t/+ ^1) *l/~"+P  [  (n  +  v  +  p  +  l)(t/  +  p  -  n) 

x  (n  +p  —  i/)(n  +  v-  p+  1)]1/J 


(A.6) 


(A.7) 
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and 


I  jl  32  J3 

y  mi  m2  -(mi  -f  m2) 

are  the  Wigner  3j  symbols,  hp(kr)  is  spherical  Hankel  function  of  first  kind  of  order 
p,  and  Y™(d,<p),  Y™(8,  <f>)  =  P™(cos8)etm*,  is  spherical  harmonic  and  P™(cos9),  m  = 
0,  ±1,  •  •  • ,  ±n  are  associated  Legendre  polynomials.  Also 


(A.8) 


(2  n  +■  l)(n  —  m)! 
47rn(n  -f  l)(n  +  m)! 


(A.9) 
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Appendix  B 


The  vector  spherical  waves  are  [12] 


RgMmn(kr)  =  7m  njn(kr)Umn($,<l>) 


(B.l) 


-  (B.2) 


where 


=  [^— j*—  ^^(-=°s»)l  »im* 

da  sin  a 


Cmn(6 


,*)  =  [*4^ 

sin0 


/>nm(cos0)-<£ 


dPnm(cosfl)]  m. 


(B.3) 

(B.4) 


(B.5) 


are  vector  spherical  harmonics  and  jfn  is  spherical  Bessel  function  of  order  n.  In  (Bl)  and 
(B2),  Rg  stands  for  regular  part  of  the  wave  function.  For  vector  waves  Mmn(kf )  and 
Nmn(kr)  without  “Rgn,  they  are  the  same  expressions  as  (Bl)  and  (B2)  respectively  with 
jn  replaced  by  hn,  spherical  Hankel  function  of  first  kind  and  of  order  n. 

The  T  matrix  coefficients  for  a  sphere  with  radius  a,  permittivity  e,,  and  wavenumber 
k,=(jjy/Jira  are 

rp(M)  _  ~~  jn(ka)[kaajn{kta)\ 

in  —  *  / 1  -WL-L  /  U  f  L  —  \IL  fL  -11/ 


n  Jn(k,a)[kahn(ko)]/  -  hn(ka)[ktajn(kta)\' 

t(N)  [kWMk'a)][kaj"(ka)]'  -  [fcaaaJn(fco)][fc.QJn(fc«a)]/ 

n  [k2a2jn{k,a)][kahn(ka)]'  -  [k2a2hn{ka)\[k.ajn{kta)]' 


(B.7) 
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Figure  I.  Examples  of  3-dimensional  aggregates  obtained  by  the  hierarchical  CCA  model, 
(a)  16  particles;  (b)  32  particles;  (c)  64  particles;  (d)  128  particles;  (e)  256  particles;  and 
(f)  512  particles. 


Figure  3(a).  Bistatic  scattering  cross  section  as  a  function  of  <f>t  for  an  aggregate  with  32 
spheres  (see  Figure  1(b)).  Other  parameters  are  /t,  =  1.75-f-t0.5,  ka  =  0.05,  0,  =  9O°,  <j>i  —  0" 
and  9,  —  90°. 


(  >  ki ) 


<t>9 


Figure  3(b).  Bistatic  scattering  cross  section  as  a  function  of  <f>t  for  an  aggregate  with  32 
spheres  (see  Figure  1(b)).  Other  parameters  are  fi,  =  1.75+t'0.5,  ka  =  0. 25,  0i  =  9O°,  <j>i=  0° 
and  9,  =  90°. 


I 
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ka  =  0.05 
fj-t  =  1.75  +  z0.5 
9i  =  90°,  &  =  O' 


Figure  4(a).  Bistatic  scattering  cross  section  as  a  function  of  <f>t  for  an  aggregate  with  128 
spheres  (see  Figure  1(d))*  Other  parameters  are  fi,  =  1.75+if  *,  Asa  =  0.05,  0,  =  9O®,  <f>i  =  0 
and  6,  =  90®. 


(Ajs,  A^) 


Figure  4(b).  Bistatic  scattering  cross  section  as  a  function  of  <pt  for  an  aggregate  with  128 
spheres  (see  Figure  1(d)).  Other  parameters  are  p,  =  1.75+*0.5,  &o  =  0.25,  0,  =  9 0",  4>i  =  0° 
and  9,  —  90°. 


1  .20 


Figure  6.  Absorption,  cross  section  as  a  function  of  number  of  primary  particles  within  a 
cluster.  Other  parameters  are  /i,  =  1.75+*0.5,  fca  =  0.05,and0.25,  0<  =  9OO,  and  <f>i  =  0°. 
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Figure  7.  Extinction  cross  section  as  a  function  of  number  of  primary  particles  within  a 
cluster.  Other  parameters  are  p,  =  1.75+i0.5,  Jbo=0.05,and0.25,  =  90°,  and  &=  0°. 


Extinction  Cross  Section  /  N 


10-3 


*  N  =  51  2 


O  =  0.03  micron 


10-7 


1  0-8 


10°  1 01  1 02 
WAVELENGTH  microns 


Figure  8.  Extinction  cross  section  per  constituent  particle  as  a  function  of  wavelength. 
The  number  of  particles  in  an  aggregate  are  N  =  128, 256,  and512.  The  radius  of  sphere  is 
a=0.03jim. 


r/(2o) 

Figure  9.  Percus-Yevick  Pair  distribution  function  for  sticky  spheres  with  r  =  0.2  and 
r  s=  0.5,  and  for  non-sticky  spheres. 
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Figure  10.  Percus-Yevick  Pair  distribution  function  for  sticky  spheres  with  /  =  0.2  and 
/  =  0.4  and  stickiness  r  =  0.2. 
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Figure  11.  Pair  distribution  function  for  media  with  sticky  spheres.  The  simulation 
parameters  are  /  =  0.2,  t  =  0.2,  N  =  512,  and  30  realizations.  The  simulated  results  are 
compared  with  Percus-Yevick  calculations. 
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Figure  13.  Radius  of  gyration  Rg  as  a  function  of  the  number  of  spherules  N.  The  squares 
are  results  of  the  SHS  model. 
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Figure  14.  Attenuation  rates  based  on  QCA-CP  and  independent  scattering  as  a  function 
of  fractional  volume  for  various  stickiness  parameters.  Other  parameters  are  t,  —  3.2ea  and 
fca  =  0.314. 
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Figure  15.  Turbidity  as  a  function  of  fractional  volume  calculated  with  the  parameters  in 
Table  2  using  QCA-SHS. 
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